################################################################################
# Immigrants, elections, and turnout
# Author: Linuz Aggeborn, Henrik Andersson, Sirus Dehdari, Karl-Oskar Lindgren
# Description: Create figures for the banwidth plots using "honest" confidence
# intervals, based on Koles�r and Rothe (2018) and Armstrong and Koles�r (2020)
################################################################################

# We save paths for input and output:
in_data <- "E://ProjData//Concurrent elections//"

out_figure <- "C://Userdata//Shared//Output//Concurrent elections//img//" 
out_table <- "C://Userdata//Shared//Output//Concurrent elections//tex//" 
out_values <- "C://Userdata//Shared//Output//Concurrent elections//values//" 

out_figure_mydp <- "//tsclient//C//Users//HP ZBook//Dropbox//Projects//Medborgarskapsprojektet//figures//"


library(readstata13) # To read Stata13 dta-files
library(xtable)
library(sandwich)
library(lmtest)
library(ggplot2)
library(Hmisc)
library(RDHonest)
library(data.table)
library(haven)



### Ordinary rounding ###
# This function will round up the values without removing last zeros (after decimal)
printround <- function(tval, round = 3, lessthanzero = FALSE){
  spr <- paste(c("%.",as.character(round),"f"), collapse="")
  tval <- sprintf(spr,round(tval,round))
  if (lessthanzero == TRUE){
    if (as.numeric(tval) == 0) tval <- paste("<0.",paste(rep("0",round-1), collapse = ""),"1", sep = "")
  }
  return(tval)
}



### Importing data: ----
finaldata <-  read.dta13(paste(in_data,"finaldata//finaldata_immigrants.dta", sep =""))


#### Create separate data for 1994 and 2010:
finaldata2010 <- finaldata[finaldata$year_ctz == 2010,]
# Create the donut:
finaldata2010 <- finaldata2010[finaldata2010$days_treat_enrolled > 7 | finaldata2010$days_treat_enrolled <=0,]
finaldata2010 <- finaldata2010[finaldata2010$days_treat_enrolled >=-30,]

finaldata1994 <- finaldata[finaldata$year_ctz == 1994,]
# Create the donut:
finaldata1994 <- finaldata1994[finaldata1994$days_treat_enrolled > 7 | finaldata1994$days_treat_enrolled <=0,]

# Creating matrices for saving estimates for 1994 and 2010, sharp RD:
bw_v <- 11:79
res94 <- matrix(ncol = 3, nrow = length(bw_v))
res94 <- cbind(bw_v, res94)
res10 <- matrix(ncol = 3, nrow = length(bw_v))
res10 <- cbind(bw_v, res10)

# Estimate RD coefficient using 11 to 79 days bandwidth:
for (i in 1:length(bw_v)){
  bw <- bw_v[i]
  tempRD94 <- RDHonest(voted_local ~ days_treat_enrolled, data = finaldata1994, h = bw, M = 0, kern = "uniform")
  res94[i,2] <- tempRD94$coefficients$estimate
  res94[i,3] <- tempRD94$coefficients$conf.low
  res94[i,4] <- tempRD94$coefficients$conf.high
  if (bw < 30){
    tempdata2010 <- finaldata2010[finaldata2010$days_treat_enrolled <=bw,]
    tempRD10 <- RDHonest(voted_local ~ days_treat_enrolled, data = tempdata2010, h = 30, M = 0, kern = "uniform")
    res10[i,2] <- tempRD10$coefficients$estimate
    res10[i,3] <- tempRD10$coefficients$conf.low
    res10[i,4] <- tempRD10$coefficients$conf.high
    }
    
  if (bw>29){  
    tempRD10 <- RDHonest(voted_local ~ days_treat_enrolled, data = finaldata2010, h = bw, M = 0, kern = "uniform")
    res10[i,2] <- tempRD10$coefficients$estimate
    res10[i,3] <- tempRD10$coefficients$conf.low
    res10[i,4] <- tempRD10$coefficients$conf.high
  }
}

coln <- c("Bandwidth", "Coefficient", "Lower", "Upper")
colnames(res94) <- coln
colnames(res10) <- coln
res94 <- data.frame(res94)
res10 <- data.frame(res10)

# Saving optimal bandwidth to be displayed in plots:
opt94 <- 44
opt10 <- 31

### Shart RDD: ----

### Bandwidth plot for 1994:
bandwidth_plot <- ggplot(res94, aes(x = Bandwidth))

bandwidth_plot + theme_bw(base_size = 20) +
  geom_line(aes(y = Coefficient), size = 1) +
  geom_line(aes(y = Lower), size = 1, linetype = "dashed") +
  geom_line(aes(y = Upper), size = 1, linetype = "dashed") +
  # geom_vline(xintercept = opt94, linetype = "longdash", color = "black") + 
#  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  xlab("Symmetric bandwidth in days") + ylab("Estimated coefficient") 
# coord_cartesian(xlim = c(bwidth_vec[1], bwidth_vec[length(bwidth_vec)]), expand = FALSE)

ggsave(paste(out_figure,"bandbreddgraf_honestCI_1994.pdf", sep = ""))
ggsave(paste(out_figure_mydp,"bandbreddgraf_honestCI_1994.pdf", sep = ""))


### Bandwidth plot for 2010:
bandwidth_plot <- ggplot(res10, aes(x = Bandwidth))

bandwidth_plot + theme_bw(base_size = 20) +
  geom_line(aes(y = Coefficient), size = 1) +
  geom_line(aes(y = Lower), size = 1, linetype = "dashed") +
  geom_line(aes(y = Upper), size = 1, linetype = "dashed") +
  # geom_vline(xintercept = opt10, linetype = "longdash", color = "black") + 
  #  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  xlab("Asymmetric bandwidth in days, upper limit") + ylab("Estimated coefficient") 
# coord_cartesian(xlim = c(bwidth_vec[1], bwidth_vec[length(bwidth_vec)]), expand = FALSE)

ggsave(paste(out_figure,"bandbreddgraf_honestCI_2010.pdf", sep = ""))
ggsave(paste(out_figure_mydp,"bandbreddgraf_honestCI_2010.pdf", sep = ""))



# Creating matrices for the 1994 and 2010 fuzzy RD estimates:
fzres94 <- matrix(ncol = 3, nrow = length(bw_v))
fzres94 <- cbind(bw_v, fzres94)

fzres10 <- matrix(ncol = 3, nrow = length(bw_v))
fzres10 <- cbind(bw_v, fzres10)

# Estimate fuzzy RD coefficient using 10 to 79 days bandwidth:
for (i in 1:length(bw_v)){
  bw <- bw_v[i]
  tempRD94 <- RDHonest(voted_local ~ eligeble | days_treat_enrolled, data = finaldata1994, h = bw, M = c(0,0), kern = "uniform")
  fzres94[i,2] <- tempRD94$coefficients$estimate
  fzres94[i,3] <- tempRD94$coefficients$conf.low
  fzres94[i,4] <- tempRD94$coefficients$conf.high
  if (bw < 30){
    tempdata2010 <- finaldata2010[finaldata2010$days_treat_enrolled <=bw,]
    tempRD10 <- RDHonest(voted_local ~ eligeble | days_treat_enrolled, data = tempdata2010, h = 30, M = c(0,0), kern = "uniform")
    fzres10[i,2] <- tempRD10$coefficients$estimate
    fzres10[i,3] <- tempRD10$coefficients$conf.low
    fzres10[i,4] <- tempRD10$coefficients$conf.high
  }
  if (bw>29){  
    tempRD10 <- RDHonest(voted_local ~ eligeble | days_treat_enrolled, data = finaldata2010, h = bw, M = c(0,0), kern = "uniform")
    fzres10[i,2] <- tempRD10$coefficients$estimate
    fzres10[i,3] <- tempRD10$coefficients$conf.low
    fzres10[i,4] <- tempRD10$coefficients$conf.high
  }
}


coln <- c("Bandwidth", "Coefficient", "Lower", "Upper")
colnames(fzres94) <- coln
colnames(fzres10) <- coln

fzres94 <- data.frame(fzres94)
fzres10 <- data.frame(fzres10)


### Fuzzy RDD: ----

### Bandwidth plot for 1994, fuzzy:
bandwidth_plot <- ggplot(fzres94, aes(x = Bandwidth))

bandwidth_plot + theme_bw(base_size = 20) +
  geom_line(aes(y = Coefficient), size = 1) +
  geom_line(aes(y = Lower), size = 1, linetype = "dashed") +
  geom_line(aes(y = Upper), size = 1, linetype = "dashed") +
  # geom_vline(xintercept = opt94, linetype = "longdash", color = "black") + 
  #  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  xlab("Symmetric bandwidth in days") + ylab("Estimated coefficient") 
# coord_cartesian(xlim = c(bwidth_vec[1], bwidth_vec[length(bwidth_vec)]), expand = FALSE)

ggsave(paste(out_figure,"bandbreddgraf_honestCI_1994_fuzzy.pdf", sep = ""))
ggsave(paste(out_figure_mydp,"bandbreddgraf_honestCI_1994_fuzzy.pdf", sep = ""))


### Bandwidth plot for 2010, fuzzy:
bandwidth_plot <- ggplot(fzres10, aes(x = Bandwidth))

bandwidth_plot + theme_bw(base_size = 20) +
  geom_line(aes(y = Coefficient), size = 1) +
  geom_line(aes(y = Lower), size = 1, linetype = "dashed") +
  geom_line(aes(y = Upper), size = 1, linetype = "dashed") +
  # geom_vline(xintercept = opt10, linetype = "longdash", color = "black") + 
  #  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  xlab("Asymmetric bandwidth in days, upper limit") + ylab("Estimated coefficient") 
# coord_cartesian(xlim = c(bwidth_vec[1], bwidth_vec[length(bwidth_vec)]), expand = FALSE)

ggsave(paste(out_figure,"bandbreddgraf_honestCI_2010_fuzzy.pdf", sep = ""))
ggsave(paste(out_figure_mydp,"bandbreddgraf_honestCI_2010_fuzzy.pdf", sep = ""))
















